## Flynn Stewart study 1 (Qualtrics, US sample, Feb. 2016)

#setwd("~/Downloads/") #set your working directory here
library(foreign)
data<-read.csv("FlynnStewart study 1 data.csv")
names(data)
dim(data) #n=605

## Recodes:
library(car)

###########################
### Treatment indicator ###
###########################

data$Q188<-recode(data$Q188,"1=1;else=0",as.factor.result=F)
table(data$Q188)
data$Q257<-recode(data$Q257,"1=1;else=0",as.factor.result=F)
table(data$Q257)
data$Q258<-recode(data$Q258,"1=1;else=0",as.factor.result=F)
table(data$Q258)
data$Q259<-recode(data$Q259,"1=1;else=0",as.factor.result=F)
table(data$Q259)
data$Q260<-recode(data$Q260,"1=1;else=0",as.factor.result=F)
table(data$Q260)
data$Q350<-recode(data$Q350,"1=1;else=0",as.factor.result=F)
table(data$Q350)

data$condition<-ifelse(data$Q188==1, 1,
                       ifelse(data$Q257==1, 2,
                              ifelse(data$Q258==1, 3,
                                     ifelse(data$Q259==1, 4,
                                            ifelse(data$Q260==1, 5,
                                                   ifelse(data$Q350==1, 6, NA))))))

table(data$condition) #Ns: 107, 114, 84, 99, 105, 96

# Approval of tactics DV
data$approve.tactics<-recode(data$Q342,"1=1;2=3;3=4;4=5;5=6;6=2;8=7;else=NA",as.factor.result=F)
table(data$approve.tactics)
tapply(data$approve.tactics, data$condition, mean, na.rm=T)

t.test(data$approve.tactics[data$condition==1, 
       data$approve.tactics[data$condition==2], na.rm=T) #

t.test(data$approve.tactics[data$condition==3, 
       data$approve.tactics[data$condition==4], na.rm=T) #

t.test(data$approve.tactics[data$condition==5, 
       data$approve.tactics[data$condition==6], na.rm=T) #

# Overall group approval DV
data$approve.group<-recode(data$Q343,"1=1;2=2;3=4;=5;10=3;13=6;14=7;else=NA",as.factor.result=F)
table(data$approve.group)
tapply(data$approve.group, data$condition, mean, na.rm=T)

t.test(data$approve.group[data$condition==1], 
       data$approve.group[data$condition==2], na.rm=T) #
       
       t.test(data$approve.group[data$condition==3], 
       data$approve.group[data$condition==4], na.rm=T) #
       
       t.test(data$approve.group[data$condition==5], 
       data$approve.group[data$condition==6], na.rm=T) #

# KNU feeling thermometer DV
data$ft<-recode(data$Q347,"NA=NA",as.factor.result=F)
summary(data$ft)
hist(data$ft)
tapply(data$ft, data$condition, mean, na.rm=T)

t.test(data$ft[data$condition==1], data$ft[data$condition==2], na.rm=T)
t.test(data$ft[data$condition==3], data$ft[data$condition==4], na.rm=T)
t.test(data$ft[data$condition==5], data$ft[data$condition==6], na.rm=T)

# KNU a legitimate alternative to Myanmar gov't? DV
data$legitimate<-recode(data$Q346, "1=1;2=2;3=3;4=4;5=5;6=6;7=7;else=NA")
table(data$legitimate)
tapply(data$legitimate, data$condition, mean, na.rm=T)

t.test(data$legitimate[data$condition==1], data$legitimate[data$condition==2], na.rm=T) #p<.01
t.test(data$legitimate[data$condition==3], data$legitimate[data$condition==4], na.rm=T) #p=.21
t.test(data$legitimate[data$condition==5], data$legitimate[data$condition==6], na.rm=T) #p=.86

# KNU capable of governing independent state? DV
data$capable<-recode(data$Q344, "2=1;5=2;6=3;7=4;8=5;else=NA",as.factor.result=F)
table(data$capable)
tapply(data$capable, data$condition, mean, na.rm=T)

t.test(data$capable[data$condition==1], data$capable[data$condition==2], na.rm=T) #p<.01
t.test(data$capable[data$condition==3], data$capable[data$condition==4], na.rm=T) #p=.86
t.test(data$capable[data$condition==5], data$capable[data$condition==6], na.rm=T) #p=.10

# Support independent KNU state? DV
data$independent<-recode(data$Q345,"1=1;11=2;12=3;13=4;14=5;15=6;16=7;else=NA",as.factor.result=F)
table(data$independent)
tapply(data$independent, data$condition, mean, na.rm=T)

t.test(data$independent[data$condition==1], data$independent[data$condition==2], na.rm=T) #p=.03
t.test(data$independent[data$condition==3], data$independent[data$condition==4], na.rm=T) #p=.13
t.test(data$independent[data$condition==5], data$independent[data$condition==6], na.rm=T) #p=.13


#########################################
### Covariates / randomization checks ###
#########################################

## Demographics/political:
data$age<-(2016 - data$Q181)
table(data$age)
data$education<-recode(data$Q182,"1=1;2=2;3=3;4=4;5=5;else=NA",as.factor.result=F)
table(data$education)
data$male<-recode(data$Q15,"1=1;else=0",as.factor.result=F)
table(data$male)
data$female<-ifelse(data$male==1,0,1)
table(data$female)
data$white.temp<-unclass(data$dem_racecps_white)
table(data$white.temp)
data$white<-recode(data$Q16,"1=1;else=0",as.factor.result=F)
table(data$white)
#PID:
data$pid.7<-recode(data$Q90,"1=1;2=2;3=3;4=4;5=5;6=6;7=7;else=NA",as.factor.result=F)
table(data$pid.7)
data$dem<-recode(data$pid.7,"1:3=1;4:7=0;else=NA",as.factor.result=F)
table(data$dem)
data$rep<-recode(data$pid.7,"1:4=0;5:7=1;else=NA",as.factor.result=F)
table(data$rep)
data$ind<-recode(data$pid.7,"1:3=0;4=1;5:7=0;else=NA",as.factor.result=F)
table(data$ind)
data$ideo.7<-recode(data$Q37,"1=1;2=2;3=3;4=4;5=5;6=6;7=7;else=NA",as.factor.result=F)
table(data$ideo.7)
data$interest<-recode(data$Q54,"1=5;2=4;3=3;4=2;5=1;else=NA",as.factor.result=F)
table(data$interest) 
data$days_talk.temp<-unclass(data$discuss_discpstwk)
data$days_talk<-recode(data$days_talk.temp,"6=0;7=1;8=2;9=3;10=4;11=5;12=6;13=7;else=NA",as.factor.result=F)
table(data$days_talk)

## General pol. knowledge (count skips as incorrect):
data$pres_terms<-recode(data$Q38,"2=1;'2 consecutive terms'=1;
                        '2 terms'=1;'2 terms 8 yrs'=1;'2 times'=1;
                        '8 years'=1;'twice'=1;'Twice'=1;'Twice for a total of eight years'=1;
                        'Twice two 4 year terms'=1;'Twice. 2 terms'=1;'two'=1;
                        'two '=1;'twoce'=1;else=0",as.factor.result=F)
table(data$pres_terms)
data$deficit<-recode(data$Q39,"1=1;else=0",as.factor.result=F)
table(data$deficit)
data$sen_term<-recode(data$Q40,"6=1;'6 years'=1;'6 years '=1;'6yrs'=1;
                      'six'=1;'Six'=1;else=0",as.factor.result=F)
data$medicare<-recode(data$Q41,"1=1;else=0",as.factor.result=F)
table(data$medicare)
data$spend_least<-recode(data$Q42,"1=1;else=0",as.factor.result=F)
table(data$spend_least)
data$speaker<-recode(data$Q185,"1=1;else=0",as.factor.result=F)
table(data$speaker)
data$vp<-recode(data$Q44,"2=1;else=0",as.factor.result=F)
table(data$vp)

data$know.score <- (data$pres_terms+
                      data$deficit+
                      data$sen_term+
                      data$medicare+
                      data$spend_least+
                      data$speaker+
                      data$vp)
summary(data$know.score)
hist(data$know.score)
median(data$know.score, na.rm=T)

#define high knowledge as above median (5 correct out of 9):
data$high.know<-ifelse(data$know.score>median(data$know.score,na.rm=T),1,0)
table(data$high.know)

## Rand checks:
tapply(data$age, data$condition, median, na.rm=T)
chisq.test(data$age,data$condition) #p=.19

tapply(data$education, data$condition, median, na.rm=T)
chisq.test(data$education,data$condition) #p=.40

tapply(data$male, data$condition, mean, na.rm=T)
chisq.test(data$male,data$condition) #p=.07

tapply(data$white, data$condition, mean, na.rm=T)
chisq.test(data$white,data$condition) #p=.60

tapply(data$dem, data$condition, mean, na.rm=T)
chisq.test(data$dem,data$condition) #p=.48

tapply(data$rep, data$condition, mean, na.rm=T)
chisq.test(data$rep,data$condition) #p=.68

tapply(data$ideo.7, data$condition, mean, na.rm=T)
chisq.test(data$ideo.7,data$condition) #p=.44

## Regressions for paper (Table 1): 
#create violence and service dummies:
data$violent<-ifelse(data$condition==2 | data$condition==4 | data$condition==6, 1, 0)
table(data$violent, data$condition)
data$restrictive<-ifelse(data$condition==3 | data$condition==4, 1, 0)
table(data$restrictive, data$condition)
data$inclusive<-ifelse(data$condition==5 | data$condition==6, 1, 0)
table(data$inclusive, data$condition)

reg1<-lm(legitimate~violent*restrictive + violent*inclusive, data=data)
summary(reg1)
reg2<-lm(capable~violent*restrictive + violent*inclusive, data=data)
summary(reg2)
reg3<-lm(independent~violent*restrictive + violent*inclusive, data=data)
summary(reg3)

##testing restrictive vs. inclusive in all three models:
linearHypothesis(reg1, "restrictive = inclusive") #n.s.
linearHypothesis(reg2, "restrictive = inclusive") #n.s.
linearHypothesis(reg3, "restrictive = inclusive") #n.s.

##testing treatment effects on dichotomized outcomes (legitimate/not legitimate:)
data$legitimate.b<-ifelse(data$legitimate %in% c(5,6,7),1,0)
table(data$legitimate,data$legitimate.b)
data$capable.b<-ifelse(data$capable %in% c(4,5),1,0)
table(data$capable,data$capable.b)
data$independent.b<-ifelse(data$independent %in% c(5,6,7),1,0)
table(data$independent,data$independent.b)

reg4<-lm(legitimate.b~violent*restrictive + violent*inclusive, data=data)
summary(reg4)
reg5<-lm(capable.b~violent*restrictive + violent*inclusive, data=data)
summary(reg5)
reg6<-lm(independent.b~violent*restrictive + violent*inclusive, data=data)
summary(reg6)

reg7<-glm(legitimate.b~violent*restrictive + violent*inclusive, data=data, family=binomial(link="logit"))
summary(reg7)
reg8<-glm(capable.b~violent*restrictive + violent*inclusive, data=data, family=binomial(link="logit"))
summary(reg8)
reg9<-glm(independent.b~violent*restrictive + violent*inclusive, data=data, family=binomial(link="logit"))
summary(reg9)

##Confidence intervals for Violence x Services interaction effects:
#NOTE: confidence intervals around these effects were used to make error bars in Figure 1
#make two treatments factors:
data$violence.fac<-ifelse(data$condition==1 | data$condition==3 | data$condition==5, "base-nv", "v")
data$violence.fac<-as.factor(data$violence.fac)
data$service.fac<-ifelse(data$condition==1 | data$condition==2, "base-none",
                         ifelse(data$condition==3 | data$condition==4, "rest", 
                                ifelse(data$condition==5 | data$condition==6, "incl", NA)))
data$service.fac<-as.factor(data$service.fac)
table(data$condition, data$violence.fac) #confirm 
table(data$condition, data$service.fac) #confirm
#Marginal effects of (violence | services) on LEGITIMACY DV:
reg1<-lm(legitimate~violence.fac*service.fac, data=data)
summary(reg1) #coefs for Table 3 in paper
library(contrast) #calculate marginal effects using "contrast" package:
#Effect of violence among groups that provide NO services:
ViolenceEffect_NoServices <- contrast(reg1,
list(service.fac = "base-none", violence.fac = "v"),
list(service.fac = "base-none", violence.fac = "base-nv"))
print(ViolenceEffect_NoServices, X = T) #Contrast=effect of violence 
#Effect of violence among groups that provide RESTRICTIVE services:
ViolenceEffect_RestServices <- contrast(reg1,
                                      list(service.fac = "rest", violence.fac = "v"),
                                      list(service.fac = "rest", violence.fac = "base-nv"))
print(ViolenceEffect_RestServices, X = T) #Contrast=effect of violence 
#Effect of violence among groups that provide INCLUSIVE services:
ViolenceEffect_InclServices <- contrast(reg1,
                                        list(service.fac = "incl", violence.fac = "v"),
                                        list(service.fac = "incl", violence.fac = "base-nv"))
print(ViolenceEffect_InclServices, X = T) #Contrast=effect of violence

#looking at extent of attenuation in violence effect across service provision 
#(t-test: violent no services vs. violent inclusive services):
t.test(data$legitimate[data$violent==1 & data$restrictive==0 & data$inclusive==0],
       data$legitimate[data$violent==1 & data$restrictive==0 & data$inclusive==1], na.rm=T)



#Marginal effects of (violence | services) on CAPABLE OF GOVERNING DV:
reg2<-lm(capable~violence.fac*service.fac, data=data)
summary(reg2) #coefs for Table 3 in paper
#Effect of violence among groups that provide NO services:
ViolenceEffect_NoServices <- contrast(reg2,
                                      list(service.fac = "base-none", violence.fac = "v"),
                                      list(service.fac = "base-none", violence.fac = "base-nv"))
print(ViolenceEffect_NoServices, X = T) #Contrast=effect of violence 
#Effect of violence among groups that provide RESTRICTIVE services:
ViolenceEffect_RestServices <- contrast(reg2,
                                        list(service.fac = "rest", violence.fac = "v"),
                                        list(service.fac = "rest", violence.fac = "base-nv"))
print(ViolenceEffect_RestServices, X = T) #Contrast=effect of violence 
#Effect of violence among groups that provide INCLUSIVE services:
ViolenceEffect_InclServices <- contrast(reg2,
                                        list(service.fac = "incl", violence.fac = "v"),
                                        list(service.fac = "incl", violence.fac = "base-nv"))
print(ViolenceEffect_InclServices, X = T) #Contrast=effect of violence

#Marginal effects of (violence | services) on SUPPORT INDEPENDENT STATE DV:
reg3<-lm(independent~violence.fac*service.fac, data=data)
summary(reg3) #coefs for Table 3 in paper
#Effect of violence among groups that provide NO services:
ViolenceEffect_NoServices <- contrast(reg3,
                                      list(service.fac = "base-none", violence.fac = "v"),
                                      list(service.fac = "base-none", violence.fac = "base-nv"))
print(ViolenceEffect_NoServices, X = T) #Contrast=effect of violence 
#Effect of violence among groups that provide RESTRICTIVE services:
ViolenceEffect_RestServices <- contrast(reg3,
                                        list(service.fac = "rest", violence.fac = "v"),
                                        list(service.fac = "rest", violence.fac = "base-nv"))
print(ViolenceEffect_RestServices, X = T) #Contrast=effect of violence 
#Effect of violence among groups that provide INCLUSIVE services:
ViolenceEffect_InclServices <- contrast(reg3,
                                        list(service.fac = "incl", violence.fac = "v"),
                                        list(service.fac = "incl", violence.fac = "base-nv"))
print(ViolenceEffect_InclServices, X = T) #Contrast=effect of violence


##Export regression tables:
library(stargazer)
stargazer(reg1, reg2, reg3, type="latex",  title="OLS Regression Models Predicting Legitimacy Evaluations (Study 1)",
          covariate.labels=c("Violent", "Restrictive Services","Inclusive Services", "Violent x Restrictive Services","Violent x Inclusive Services","Constant"), omit.stat=c("rsq", "adj.rsq", "res.dev", "F"),
          dep.var.labels = c("Legitimate Alternative","Capable of Governing","Support Independent State"),
          model.numbers=F, digits=2)

#robustness check: dichotomized outcomes with OLS (for appendix):
stargazer(reg1, reg2, reg3, reg4, reg5, reg6, type="latex",  title="Regression Models Predicting Legitimacy Evaluations (Study 1)",
          covariate.labels=c("Violent", "Restrictive Services","Inclusive Services", "Violent x Restrictive Services","Violent x Inclusive Services","Constant"), omit.stat=c("rsq", "adj.rsq", "res.dev", "F"),
          dep.var.labels = c("Legitimate Alternative (1-7)","Capable of Governing (1-5)","Support Independent State (1-7)",
                             "Legitimate Alternative (0-1)", "Capable of Governing (0-1)", "Support Independent State (0-1)"),
          model.numbers=F, digits=2)

#robustness check: dichotomized outcomes with logit (for appendix):
stargazer(reg1, reg2, reg3, reg7, reg8, reg9, type="latex",  title="Regression Models Predicting Legitimacy Evaluations (Study 1)",
          covariate.labels=c("Violent", "Restrictive Services","Inclusive Services", "Violent x Restrictive Services","Violent x Inclusive Services","Constant"), omit.stat=c("rsq", "adj.rsq", "res.dev", "F"),
          dep.var.labels = c("Legitimate Alternative (1-7)","Capable of Governing (1-5)","Support Independent State (1-7)",
                             "Legitimate Alternative (0-1)", "Capable of Governing (0-1)", "Support Independent State (0-1)"),
          model.numbers=F, digits=2)

#Export additional table with controls for appendix (model 1 controlling for gender; model 2 controlling for age, education, gender, race, Dem, Rep, ideology):
reg10<-lm(legitimate~violent+restrictive+inclusive+violent:restrictive+violent:inclusive+male, data=data)
summary(reg10)
reg11<-lm(legitimate~violent+restrictive+inclusive+violent:restrictive+violent:inclusive+male+age+education+white+dem+rep+ideo.7, data=data)
summary(reg11)

stargazer(reg8, reg9, type="latex",  title="Robustness Checks (Study 1)",
          dep.var.labels = c("Legitimate Alternative"),
          model.numbers=F, digits=2) #note: changed order of coefs manually in Latex
